Biogeochemical niche conservatism relates to plant species diversification and life form evolution in a subtropical montane evergreen broad‐leaved forest

Abstract The evolutionary mechanisms underlying the biogeochemical niche conservatism in forests remain incompletely understood. Here we aimed to determine how the strengths of biogeochemical niche conservatism vary among elements and between life forms. We measured leaf concentrations of basal elements (C, N, P, K, Ca, and Mg) in a wide range of life forms in a subtropical montane evergreen broad‐leaved forest. We found that differences in life forms such as evergreen/deciduous woody species and herbaceous/woody species significantly affected leaf elemental composition. The significant phylogenetic signal was present in leaf C, N, K, and Mg concentrations but absent in leaf P and Ca concentrations in all species. These contrasting strengths of biogeochemical niche conservatism were best generated by Ornstein–Uhlenbeck processes toward optima. Woody species were evolutionarily selected to show lower optimal leaf N, P, and K concentrations and higher optimal leaf C, Ca, and Mg concentrations than herbaceous species. The number of optima varied from the least in leaf C concentration to the most in leaf Ca concentration, suggesting the stronger convergent evolution of leaf Ca concentration. The positions of optima toward the tips were more selected in woody species, suggesting the more frequency of species‐specific adaptations in woody species. The positions of optima were also selected at the nodes towards the species groupings from certain life forms (e.g., the group of 12 Polypodiales ferns in leaf Ca evolution and the group of three evergreen Theaceae species in leaf P evolution) that were converged to present similar leaf elemental composition. During the evolution of biogeochemical niche, strong correlations were found among leaf C, N, P, and K concentrations and between leaf Ca and Mg concentrations. In conclusion, the strengths of biogeochemical niche conservatism can vary among elements and between life forms due to the different tempo and mode of Ornstein–Uhlenbeck processes.


| INTRODUC TI ON
Plant species vary greatly in the quantities of essential elements in leaves to survive and grow (Ågren, 2008;de la Riva et al., 2018).
For example, the ranges of variation in leaf elemental concentrations were 2-, 32-, 75-108-, 764-, and 36-fold for C, N, P, K, Ca, and Mg in global plant species, respectively (Ma et al., 2018;Metali et al., 2015). These contrasting concentration ranges of the six essential elements in leaves across plant species suggest that species generally have different requirements of essential elements in leaves and should thus have their own leaf elemental composition, which can be called biogeochemical niche (Peñuelas et al., 2008(Peñuelas et al., , 2019. Biogeochemical niches in plant species can be phylogenetically conserved . Biogeochemical niche conservatism is usually estimated by the phylogenetic signal in leaf elemental composition, that is, the tendency for closely related species to show similar leaf elemental composition (e.g., Peñuelas et al., 2010;Sardans et al., 2021). Significant phylogenetic signal has been detected in leaf elemental composition in field-collected plants (e.g., leaf C, N, P, K, and Mg concentrations in tropical and temperate forest plant species), lending support to biogeochemical niche conservatism (Duffin et al., 2019;Peñuelas et al., 2010). However, the evolutionary mechanisms generating biogeochemical niche conservatism in plant species remain incompletely understood Peñuelas et al., 2019).
Discerning the evolutionary mechanisms responsible for biogeochemical niche conservatism in plant species is a challenge that is complicated by three problems. First, biogeochemical niche conservatism can be generated by a variety of evolutionary processes. Cooper et al. (2010) have outlined that the patterns of niche conservatism can be underlain by five evolutionary models (drift, niche retention, niche filling/shifting, phylogenetic inertia, and evolutionary rate). These models describe that species niches have evolved through a pure Brownian motion (BM) process (i.e., drift model), a more expected BM process (i.e., niche retention model), a BM process with decreased/accelerated evolutionary rate (i.e., niche filling/ shifting model), or a BM process with the incorporation of stabilizing selection (i.e., phylogenetic inertia model; Cooper et al., 2010).
Previous studies have found that the values of Blomberg's K or Pagel's λ in leaf elemental concentrations are usually found to be lower than unity (i.e., pure BM process) even if they differ significantly from zero (i.e., no phylogenetic signal; e.g., Duffin et al., 2019;Hao et al., 2015;Peñuelas et al., 2010;Sardans et al., 2015), suggesting that the generation of biogeochemical niche conservatism might not follow a drift or niche retention model (Cooper et al., 2010). Several studies have shown that the values of Pagel's δ were higher than unity for leaf N, P, K, Ca, Mg, S, Fe, and Zn concentrations in tropical karst tree species and cave-dwelling herbaceous species (Bai et al., 2019(Bai et al., , 2020, suggesting that biogeochemical niche conservatism could be generated by niche shifting model (Cooper et al., 2010;Pagel, 1999).
Alternatively, some studies have shown that the generation of biogeochemical niche conservatism might follow the phylogenetic inertia model, that is, stabilizing selection toward an optimum or optima (Bai et al., 2019(Bai et al., , 2020Fernández-Martínez et al., 2018). For example, Fernández-Martínez et al. (2018) found that the evolution of leaf P concentration could be best described by stabilizing selection towards an optimum in global plants. In contrast, they found that the higher values of optimal leaf N concentration were evolutionarily selected in some families such as Betulaceae and Salicaceae.
In addition, a high evolutionary rate of change in one group compared with another is evidence of less-constrained evolution for the same trait in two groups (i.e., evolutionary rate model ;Cooper et al., 2010). This is the case in biogeochemical niche conservatism.
For instance, Bai et al. (2020) found the lower evolutionary rates of leaf N, P, K, Ca, Mg, S, Fe, Mn, and Zn concentrations in ferns than herbs in caves. Since the generation of biogeochemical niche conservatism has been supported by different evolutionary models, researchers can compare which evolutionary model best fits dataset (Cooper et al., 2010;Münkemüller et al., 2015).
Second, the strength of biogeochemical niche conservatism might vary among elements. Previous studies have often identified that phylogenetic signals were present in some leaf elements but absent in other leaf elements in an environment. For example, Sardans et al. (2015) have detected that phylogenetic signal was significant in leaf N, K, Ca, Mg, and S concentrations but weak in leaf P concentration in European forest tree species. The contrasting levels of phylogenetic signals among leaf elements might be related to their different evolutionary rate and processes. This is supported by the finding that phylogenetic signal was strongly and negatively correlated with evolutionary rate and the strengths of niche shifting and phylogenetic inertia across leaf N, P, K, Ca, Mg, S, Fe Cu, and Zn concentrations in tropical karst tree species (Bai et al., 2019). The contrasting levels of phylogenetic signals among leaf elements might also be related to their different nutritional statuses. For example, the generally observed P limitation in tropical and subtropical soils could make distantly related species to present similar and low leaf P concentration and thus erase the phylogenetic signal in leaf P concentration (Bai et al., 2019;Metali et al., 2012;Wen et al., 2018). In addition, the evolution of leaf elements could be coordinated. A typical example is that leaf Ca and Mg concentrations are commonly observed to be positively interacted among specific groups of angiosperms when their growth is K E Y W O R D S biogeochemical niche, conservatism, convergence, element interactions, life forms, Ornstein-Uhlenbeck process

T A X O N O M Y C L A S S I F I C A T I O N
Biogeochemistry, Evolutionary ecology not limited by mineral nutrition, indicating that a significant historical relationship between the two elements exists throughout the evolutionary history of angiosperms (Hernández et al., 2013;Neugebauer et al., 2018 Bai et al. (2020) have demonstrated that herbs tend to have higher evolutionary rate than ferns for leaf N, P, K, Ca, Mg, S, Fe, Mn, and Zn concentrations in caves. Moreover, they found that the higher optimal values of leaf Ca, N, and Mg were evolutionary selected in cave-dwelling herbs. Therefore, we predicted that the contrasting strengths of biogeochemical niche conservatism between life forms might be related to their different evolutionary rates and processes in an environment.
The triptych of problems suggests the necessity to choose the best-fitted evolutionary process to evaluate the strengths of biogeochemical niche conservatism among elements and between life forms. To our knowledge, little attention has been paid to determine the evolutionary processes responsible for the contrasting biogeochemical niche conservatism among elements and between life forms. In this work, our goal was thus to determine how the strength of biogeochemical niche conservatism differs among elements and between life forms. To address this goal, we measured leaf concentrations of basal elements (C, N, P, K, Ca, and Mg) in a wide range of herbaceous and woody species in a subtropical montane evergreen broad-leaved forest in southern China (Bai et al., 2015). Similar to the typical lowland tropical forests, the soil in this subtropical forest is highly leached due to the strong mean annual precipitation (i.e., >2500 mm) and thus tends to be P-limited (Huang & Jiang, 2002;Spicer et al., 2020). However, this subtropical forest is species-rich and contains more woody species than herbaceous species (Huang & Jiang, 2002;Spicer et al., 2020). Within herbaceous species concentrated on the forest floor, ferns are more abundant than herbs. Within woody species dominating the forest understory and canopy layers, there also exists the variation in species richness due to leaf habits.
Evergreen lineages such as Fagaceae are more abundant than deciduous lineages such as Sapindaceae (Huang & Jiang, 2002 Figure A1 in Appendix). The forest is primarily old-growth but has disturbed parts due to natural disasters such as winter ice storm and human disturbances such as ecological tourism. The forest is mainly dominated by Castanopsis fabri, C. fargesii, C. caresii, Schima superba, and Rhododendron elliptic in the forest canopy and Diplopterygium chinese, Woodwardia japonica, and Lophatherum gracile in the forest floor. Soil substrate in the forest is montane yellow soil. Soil pH, concentrations of organic matter, N, P and K, and cation exchange capacity are 4.57%, 14%, 0.5%, 0.12%, 2.1% and 28.7 me/100 g, respectively (Huang & Jiang, 2002). According to the climatic observations conducted by Guangxi Lijiangyuan Forest Ecosystem Research Station near to the forest (i.e., at 1100 m altitude in Mao'er Mountain), mean annual temperature and mean annual precipitation from 2016 to 2020 were 14.0°C and 3342 mm, respectively.

| Plant sampling and measurement
We sampled the co-existing herbaceous and woody species in the forest in July 2019 and 2020. The sampled herbaceous species included 19 ferns and 16 herbs. Among woody species, the number of evergreen and deciduous species was 43 and 24, respectively.
Collectively, the total 102 plant species were from 77 genera and 45 families ( Table A1 in Appendix).
To quantify leaf elemental concentrations, we collected leaf samples using a unified protocol (Zhao et al., 2016). For herbaceous species, we sampled a certain number of individuals per species that could provide enough leaf samples for measurement.
For woody species, we sampled the current year's fully expanded and canopy leaves from three individuals per species. The harvested leaf samples of each species were oven-dried at 60°C for 72 h and ground to fine powder for chemical analyses. Each species had three measurement replicates. To determine leaf C and N concentrations, leaf samples were analyzed using an elemental analyzer (Vario MAX CN Elemental Analyzer). Leaf P, K, Ca, and Mg concentrations were determined using an inductively coupled plasma emission spectrum (iCAP Qc, Thermo Fisher Scientific; Bai et al., 2020).

| Testing for biogeochemical niche differences in plants
To test biogeochemical niche differences between life forms (i.e., ferns/herbs, evergreen/deciduous woody species, and herbaceous/ woody species), we first performed analyses of variance (ANOVA) for each element using log 10 -transformed leaf elemental concentrations.
In addition to univariate element, we also investigated the evolution of biogeochemical niche in a multivariate sense, using a principal component analysis (PCA) on mean leaf elemental concentrations across all species . The first two PCA axes explained a large proportion of total variation in leaf elemental composition. We then performed ANOVAs to detect whether the species scores of the first two PCA axes (PC1 and PC2) differed between life forms. The above conventional statistical analyses were conducted in SPSS 13.0 (SPSS).

| Testing for biogeochemical niche conservatism in plants
We constructed a phylogenetic tree of the sampled species using the V. PhyloMaker package in R (Jin & Qian, 2019) that generates phylogenies for plants based on the backbone mega-tree of 10,587 genera and 74,533 species of vascular plants, that is, GBOTB (Smith & Brown, 2018;Zanne et al., 2014). All 45 families in our study were included in GBOTB and thus were well resolved at the family level. were missing from GBOTB, they were treated as sisters and added randomly to their closely related genera in GBOTB (i.e., scenario 2 in V. PhyloMaker; Jin & Qian, 2019). We finally got a tree with distinct life forms and time-calibrated branch lengths for the subsequent phylogenetic analyses (Figure 1). Like previous studies (e.g., Zhang et al., 2021Zhang et al., , 2022, the insertion of the 18 species into the backbone tree could lead to phylogenetic uncertainty. To account for such phylogenetic uncertainty, we generated 100 random trees by scenario 2 and used them to replicate phylogenetic analyses to evaluate the stability of biogeochemical niche conservatism pattern during species diversification (Table A1 in Appendix).
To compare the strengths of biogeochemical niche conservatism among elements and between life forms, we conducted a series of phylogenetic analyses in R (version 4.1.2; R Development Core Team, 2008). Before phylogenetic analyses, species mean values of leaf elemental concentrations were log 10 -transformed and tree depth was scaled to unity via the rescale function in the GEIGER package (Harmon et al., 2008). We first detected the phylogenetic signals in leaf elemental concentrations and PC1 and PC2 species scores using Blomberg's K statistic in all species to get the general patterns of biogeochemical niche variation during species diversification (Blomberg et al., 2003). The Blomberg's K values were tested for the null hypothesis of absence of the phylogenetic signal, that is, trait values randomly distributed in the phylogenetic tree. This was achieved by randomization for K (999 times). The K values significantly larger than zero suggest significant phylogenetic signals.
The K values equal to unity suggest traits have evolved following BM process (i.e., drift model) and K values larger than unity indicate closely related species more similar than expected under BM process (i.e., niche retention model; Cooper et al., 2010). We used the phyloSignal function in the PHYLOSIGNAL package to detect phylogenetic signals (Keck et al., 2016).
We then compared a set of evolutionary models fitting to leaf elemental concentrations and PC1 and PC2 species scores in all species. Models included BM, Pagel's δ and Ornstein-Uhlenbeck (OU).
BM is a neutral model of evolution where species traits have evolved via random-walk at a constant rate (i.e., sigma). We fitted the BM models with a single evolutionary rate and two different evolutionary rates for herbaceous and woody species (BMM). We fitted BMM because the two life forms grew in contrasting light environments.
Herbaceous species were concentrated on the forest floor where low light availability might condition the evolution of biogeochemical niche. We fitted BM and BMM using the mvBM function in the mvMORPH package (Clavel et al., 2015). Pagel's δ describes that species trait evolution has accelerated or slowed down over time. The values of δ lower than unity indicate temporally early trait evolution or early-burst (i.e., niche filling model) and the values of δ larger than unity indicate temporally latter trait evolution (i.e., niche shifting model). We fitted Pagel's δ model via the fitContinuous function in the GEIGER package (Harmon et al., 2008). We tested whether δ was significantly different from unity using the pgls function in the CAPER package (Orme et al., 2012). OU describes that species traits have evolved toward an optimum or optima under a stabilizing strength (i.e., phylogenetic inertia, ɑ). We fitted the OU model with a single optimum for all species (OU1) and two optima (OUM) to test whether selective optima differed between herbaceous and woody species. We fitted OU1 and OUM using the mvOU function in the mvMORPH package. To evaluate the possible different optima within herbaceous or woody species in the phylogeny, we detected the positions of optima in the phylogeny using the OUshifts function with the singular Bayesian information criterion penalty in the PHYLOLM package (Tung Ho & Ané, 2014). Comparisons of BM, BMM, Pagel's δ, OU1, and OUM models were based on the sample size-corrected Akaike information criterion (AICc; e.g., Bai et al., 2020;Fernández-Martínez et al., 2018). When models differed by ΔAICc < 2, they were considered approximately equivalent (Burnham & Anderson, 2002), and the model with the lowest AICc value was considered the best-fitted.
To test the coordinated evolution among elements, we finally explored correlations between leaf elemental concentrations in all species by phylogenetic generalized least squares with a phylogenetic covariance matrix, simultaneously estimating Pagel's λ for phylogenetic heritability of regression residues (McCormack et al., 2020;Pagel, 1999). We used the pgls function in the CAPER package to get the element interaction strengths (Orme et al., 2012).
Leaf elemental composition significantly differed between herbaceous and woody species, which was characterized by higher leaf N, P, and K concentrations but lower leaf C, Ca, and Mg concentration F I G U R E 1 Phylogenetic relationships among the co-existing herbaceous (ferns in black, herbs in red) and woody (deciduous in blue, evergreen in green) species in a subtropical montane evergreen broad-leaved forest.

t a t h e ly p t e r is g r a c il e s c e n s C y c lo s o r u s in t e r r u p t u s H u p e r z ia s e r r a ta S e la g in e ll a u n c in a ta S e la g in e ll a ta m a ri s c in a
in herbaceous species (Table 1). Within herbaceous species, ferns tended to have lower leaf N and P concentrations than herbs (Table 1). Within woody species, leaf elemental composition was significantly affected by leaf habit (Table 1). Evergreen woody species had lower leaf N, P, K, Ca, and Mg concentrations but higher leaf C concentrations than deciduous woody species (Table 1). The PCA showed that the first two axes (i.e., PC1 and PC2) explained accumulative 73.1% of the total variation in leaf elemental composition ( Figure 2). Leaf C concentration loaded negatively and leaf N, P, and K concentrations positively with PC1 and leaf Ca and Mg concentrations loaded positively with PC2 ( Figure 2). The PC1 and PC2 species scores differed significantly between herbaceous and woody species and between evergreen and deciduous woody species (Table 1), again suggesting the biogeochemical niche variations between herbaceous and woody species and between leaf habits ( Figure 2).

| Biogeochemical niche conservatism in plants
The phylogenetic signals in leaf elemental composition (i.e., Blomberg's K) varied from 0.056 (i.e., P) to 0.245 (i.e., C) in all species (Table 2). Leaf C, N, K, and Mg concentrations had significant phylogenetic signals and leaf Ca and P concentrations exhibited non-significant phylogenetic signals. However, the phylogenetic signals in PC1 and PC2 species scores were all significant (Table 2).
Our replicated phylogenetic signal analyses on 100 random trees suggested that the frequencies of trees showing similar results to the tree in Figure 1 were 96%, 100%, 97%, 92%, 98%, 81%, 100%, 80% in leaf C, N, P, K, Ca, and Mg concentrations and PC1 and PC2 species scores, respectively (Table A1 in Appendix), indicating the highly stable pattern in biogeochemical niche conservatism during species diversification.
Among the evolutionary processes fitting to leaf elemental concentrations and PC1 and PC2 species scores in all species, BM models were least supported due to the highest AICc values ( Pagel's δ models were also better than BM models ( Table 2).
The δ values were all significantly higher than unity and varied from 46.74 in leaf C concentration to 467.08 in leaf P concentration.
OU1 models did not outperform Pagel's δ models due to their small ΔAICc differences (Table 2). But OU1 models outperformed BM models ( Table 2). The strength of phylogenetic inertia (i.e., ɑ) ranged from 23.35 in leaf C concentration to 233.03 in leaf P concentration ( Note: n is the number of species. The different upper case letters indicate significant differences between herbaceous and woody species. The different lower case letters indicate significant differences between evergreen and deciduous woody species or between ferns and herbs. leaf N, P, and K concentrations than woody species, which in turn had higher optimal values of leaf C, Ca, and Mg concentrations. Polypodiales ferns), and one node within woody species (i.e., the node toward the group of five deciduous species such as Acer dividii and Turpinia montana (Figure 3). Leaf C concentration was significantly and negatively correlated with leaf N, P, and K concentrations in all species (Table 3). Leaf Ca and Mg concentrations were significantly and positively correlated, as did for the correlations among leaf N, P, and K concentrations (Table 3).

| Evidence of biogeochemical niche differences in plants
We found large means and ranges of leaf elemental composition across all species (Table 1) can be considered as a mechanism maintaining species diversity in a given environment (Peñuelas et al., 2019;Silvertown, 2004).
We found that herbaceous and woody species had contrasting biogeochemical niches, as supported by their significantly different leaf elemental composition (Table 1). Herbaceous species in our study had significantly lower leaf C, Ca, and Mg concentrations but higher leaf N, P, and K concentrations than woody species, which is consistent with the research on savanna plants and global plants (Rossatto & Franco, 2017;Vergutz et al., 2012). This might be explained by the fact that woody species had thicker cell walls than herbaceous species, leading to higher leaf mass per area and the ensuing lower leaf N, P, and K concentrations (

F I G U R E 3
The positions (red stars) of the selective optima of species scores of the first two principal component axes (PC1 and PC2) in the phylogeny.
concentrations than deciduous woody species, suggesting the contrasting biogeochemical niches between leaf habits. This is largely due to the longer leaf life span in evergreen species that relates to larger construction cost and lower concentrations of mineral such as N, P, K, Ca, and Mg in leaves (Villar et al., 2006;Xing et al., 2021). We found that the differences of leaf elemental composition between ferns and herbs were element-specific (Table 1). Leaf C, K, Ca, and Mg concentrations differed little between ferns and herbs. However, ferns had significantly lower leaf N and P concentrations than herbs in our study, which is consistent with observations on terrestrial plant species in China (Han et al., 2005). This could be due to the strong mesophyll limitation in ferns that can reduce the requirements of N and P for photosynthesis (Tosens et al., 2016).  (Table 1), which could be a mechanism to increase carbon gain in the light-limited forest floor (Poorter et al., 2019). In contrast to herbaceous species, woody species grow in relatively light-rich environments such as sub-canopy and canopy layers where light might not be a limiting factor for plant growth. Woody species might thus develop distinct functional strategies to cope with soil and climatic factors that might constrain leaf elemental composition.
We found that leaf N, K, Ca and Mg concentrations in woody species were larger than the physiological concentration requirements (Kirkby, 2012) and bordered the sufficient concentration ranges for healthy growth (White & Brown, 2010), suggesting that N, K, Ca, and Mg might be sufficient for plant growth in our study site. However, leaf P concentration in woody species was lower than the physiological concentration requirement for P, suggesting that P might be a limiting element for plant growth. This is consistent with the general observation of P limitation in subtropical forests (Wen et al., 2018).
The P limitation, together with subtropical hot and humid climate, could facilitate evergreen woody species to develop thicker and long-lived leaves with lower mineral concentration requirements to adapt these soil and climatic conditions (Ye et al., 2022).

| Evidence of biogeochemical niche conservatism in plants
Apart from adaptations to the environmental factors, biogeochemical niches in plants can be a legacy from their ancestors, which is referred as phylogenetic conservatism (Prinzing et al., 2001). We found that leaf C, N, K, and Mg concentrations had significant phylogenetic signals in all species (  (Duffin et al., 2019). In contrast, we found that leaf Ca and P concentrations showed non-significant phylogenetic signals in subtropical forest plants ( Table 2), suggesting that phylogenetic relatedness might be less informative for the variations in the two leaf elements in subtropical forest plants. This is consistent with the detection of weak phylogenetic signal in leaf Ca and P concentrations in 58 topical lowland forest tree species (Metali et al., 2012) and in 40 tropical karst forest tree species (Bai et al., 2019). Our replicated phylogenetic signal tests on 100 randoms trees again revealed that phylogenetic conservatism was observed to be generally higher in leaf C, N, K, and Mg concentrations than leaf P and Ca concentrations (Table A1 in Appendix). We consider that these contrasting strengths of biogeochemical niche conservatism might be generated by different evolutionary processes.
The generation of biogeochemical niche conservatism was least described by BM model, because the highest AICc values were found in BM model for all leaf elements ( Table 2) between life forms. For instance, the sigma values were higher in leaf P and Ca concentrations and in woody species (Table 2), leading support to the evolutionary rate model and providing evidence of the less biogeochemical niche conservatism in the two elements and in woody species (Cooper et al., 2010;Münkemüller et al., 2015).
We found that Pagel's δ model was better than BM model and the δ values were higher than unity in all species (Table 3), which agrees with niche shifting model and indicates that the changes in biogeochemical niche in subtropical forest plants could expand to the present (Cooper et al., 2010;Pagel, 1999). The niche shifting model has also been detected to outperform BM model in the evolution of leaf elemental concentrations in tropical karst seasonal rainforest tree species (Bai et al., 2019). Under the niche shifting model, the biogeochemical niches of modern species might shift into a different portion of niche space as a consequence of selection to exploit new elemental resources, suggesting the species-specific adaptations to present environment (Cooper et al., 2010;Hernández et al., 2013).
It's noteworthy that δ value was largest in leaf P concentration, suggesting that the speed of trait trend expanding to the present was fastest in leaf P concentration. We considered that the fastest speed of trait trend expanding to the present could facilitate the speciesspecific adaptations to changes in P resources in local environments.
We found that OU1 model also outperformed BM model for all leaf elements (Table 2), supporting the phylogenetic inertia model for the evolution of biogeochemical niche in subtropical forest plants (Cooper et al., 2010). OU1 model has also been found to outperform BM model for the evolution of leaf elemental concentrations in local and global plants (e.g., Bai et al., 2019;Fernández-Martínez et al., 2018). Under OU1 model, the phylogenetic inertia (i.e., ɑ) could make leaf elemental concentrations to stabilize around the optimal composition, which could act as an evolutionary filter to exclude the species with unfit leaf elemental composition. For example, evergreen woody species with high leaf N and P concentrations might be excluded because this leaf elemental composition is not adaptive in nutrient-limited subtropical forest. We found that the ɑ value was largest in leaf P concentration (Table 3), which suggests that plant species might have the strongest phylogenetic inertia to maintain optimal leaf P concentration in subtropical forest.
We found that OUM model might be the best model to generate biogeochemical niche conservatism due to its lowest AICc values in all leaf elements (Table 2). We found that leaf elemental concentrations could evolve around the distinct optimal values between herbaceous and woody species, as shown by higher Opt H than Opt W for leaf N, P, and K concentrations ( Table 2). The higher optimal concentrations of the three elements in herbaceous species might be evolutionarily selected to widen the stomata openness (through K) for maximizing photosynthesis (through N and P), facilitating their growth and survival in light-limited forest floor (Rawat et al., 2022;Wang et al., 2013). By contrast, the higher optimal concentrations of leaf C, Ca, and Mg concentrations were selected in woody species (Table 2), leading to the greater power to maintain the stability of cell wall chemistry in leaves White et al., 2018).
We found that the number of optima under an OU process throughout the phylogeny varied from one for C to six for Ca, suggesting the different number of optima among leaf elements ( Figures A2-A7 in Appendix). According to a simulation study on trait conservatism, the increasing number of optima tends to decrease phylogenetic signal (Münkemüller et al., 2015). Consistent with this finding, leaf C concentration with the least optima presented the strongest phylogenetic signal while leaf Ca concentration with the most optima exhibited non-significant phylogenetic signal (Table 2).
This might be due to the fact that the more optima could increase the frequency of convergence, which could drive distantly related species towards the same optimal leaf elemental composition (Arbuckle et al., 2014). We found that the positions of optima could also vary among leaf elements. Some optima were concentrated towards the We found differential interaction strengths among leaf elements (Table 3). Specifically, leaf C, N and P and K concentrations were and Zn among plant species at local, regional and global scales (de la Riva et al., 2018;Watanabe et al., 2007;Zhang et al., 2012). The strong correlations among leaf elements across species might suggest that there is a phylogenetic effect governing the coordinated evolution among leaf elements (Bai et al., 2020;Donovan et al., 2011). For instance, leaf N and P concentrations have been found to be strongly and positively correlated through the phylogeny in African woody savanna species, indicating that a significant historical relationship between the two elements exists throughout the evolutionary history of this savanna plant group (Wigley et al., 2016). In addition, we found that both PC1 loaded with leaf C, N, P, and K concentrations and PC2 loaded with leaf Ca and Mg concentration presented significant phylogenetic signals (Table 2)

| CON CLUS ION
Our results demonstrate that the means and ranges of leaf elemental composition varied greatly in all species, suggestive of speciesspecific biogeochemical niche in subtropical forest plants. We found that leaf elemental composition differed significantly between life forms (i.e., herbaceous/woody species, evergreen/deciduous woody species), suggesting the important role of life form in biogeochemical niche variation. We found that phylogenetic signal was significant in leaf C, N, K, and Mg concentrations but non-significant in leaf P and Ca concentrations, which could be partly explained by the lower values of evolutionary rate (sigma), niche shifting strength (Pagel's δ) and phylogenetic inertia (α) in the former elements. However, the contrasting strengths of biogeochemical niche conservatism among elements were best generated by OU processes toward optima. We found that herbaceous species were evolutionarily selected to have higher optimal leaf N, P, K concentrations and lower optimal leaf C, Ca, and Mg concentrations than woody species, suggesting the different optima for the two life forms. We found the number of optima

ACK N OWLED G M ENT
This work was supported by National Natural Science Foundation (31960214) and Guangxi Natural Science Foundation (2020GXNSFBA297153) in China. We thank Qiumei Teng and Hao Shen for their help in plant sampling. Special thanks to the anonymous reviewers for their insightful comments on the manuscript.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data from this study are available in Table A2 in Appendix.